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ABSTRACT 

In this paper we discuss the effect of recombinations to highly excited states (« > 100) in 
hydrogen during the cosmological recombination epoch. For this purpose, we developed a 
new ODE solver for the recombination problem, based on an implicit Gear's method. This 
solver allows us to include up to 350 /-resolved shells or ~ 61 000 separate levels in the hy- 
drogen model and to solve the recombination problem for one cosmology in ~ 27 hours. This 
is a huge improvement in performance over our previous recombination code, for which a 
100-shell computation (5050 separate states) already required ~ 150 hours on a single pro- 
cessor. We show that for 350 shells down to redshift z ~ 200 the results for the free electron 
fraction have practically converged. The final modification in the free electron fraction at 
z ~ 200 decreases from about ANJN e ~ 2.8% for 100 shells to ANJN e ~ 1.6% for 350 
shells. However, the associated changes in the CMB power spectra at large multipoles / are 
rather small, so that for accurate computations in connection with the analysis of Planck data 
already ~ 100 shells are expected to be sufficient. Nevertheless, the total value of r could still 
be affected at a significant level. We also briefly investigate the effect of collisions on the re- 
combination dynamics. With our current estimates for the collisional rates we find a correction 

of AN e /N e 8.8 x 10~ 4 at z ~ 700, which is mainly caused by /-changing collisions with 

protons. Furthermore, we present results on the cosmological recombination spectrum, show- 
ing that at low frequencies collisional processes are important. However, the current accuracy 
of collisional rates is insufficient for precise computations of templates for the recombination 
spectrum at v < 1 GHz, and also the effect of collisions on the recombination dynamics suffers 
from the uncertainty in these rates. Improvements of collisional rates will therefore become 
necessary in order to obtain a final answer regarding their effects during recombination. 
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1 INTRODUCTION l200fj) . In total about 57% of all hydrogen atoms became neu- 

..... tral via the 2s- Is two-photon decay channel, while some 43% of 

Close to the maximum of the Thomson visibility function , , , , . , , T , , 

i — : — ii -i „ J the hydrogen electrons recombined through the Lyman-tt channel 

^ SunyaevfeZeldovicH B at redshift z ~ 1100, and slightly (Chk f ba & Sunvaev 2006a) . 
before, the dynamics of hydrogen recombination is mainly con- 

trolled by the net 2s-ls two-photon deca y rate and the escape At redshift z ~ 1100, the rate at which electrons meet a proton 



of photo ns fr om the Lyman-q re sonance dzeldovich et al.l [1961 : and ^combine is still very large, even though only a much smaller 



|Peebledfl96l ISunvaev & Chlubll2009ri . These two channels are fractlon of these recombinations actually end with an electron set- 

the main 'bottlenecks' during the epoch of hydrogen recombina- tlin g int0 the g round state - Under such circumstances the exact rate 

tion, with the Lyman-ff channel being more important at high red- at which electrons are captured by protons is not absolutely crucial 

shifts (z > 1300 - 1400) and the 2s- Is two -photon decay chan- for P recise computations of the cosmological ionization history, as 

nel dominating at lower redshifts (e.g. see iRubino-Martfn eTall lon 8 as this rate is close enou § h to the ^ value ' and does not 

lead to an artificial 'bottleneck' caused by the incompleteness of 

the used atomic model for hydrogen. 

* c_ tui u n •. At lower redshifts (z < 800 - 900), however, the rate of re- 

b-mail: jchluba@cita.utoronto.ca v " 

t E-mail: vasil@cita.utoronto.ca combinations drops significantly, because ever fewer free electrons 

t E-mail: ljdursi@cita.utoronto.ca and protons are available, so that the number densities of free elec- 



© 0000 RAS 



2 Chluba, Vasil <&- Dursi 



trons and protons starQ to 'freeze out'. During and just before this 
period, the total electron capture rate becomes one of the additional 
'bottlenecks' of cosmological recombination, and the precise value 
of the total recombination cross section starts to become very im- 
portant. Under these circumstances the completeness of the used 
atomic model for hydrogen becomes one of the key elements for 
accurate computations of the cosmological recombination history, 
which aim at reaching a level of precision down to ~ 0.1%. 

With measurements of the cosmic microwave background 
(CMB) temperature and polarization power spectra, as currently 
carried out with the Planck Surveyor, it indeed has become very 
important to unde rstand the ionization history at t his level of 
preci sion (e.g. see iHu et ai1ll995l : ISeliak et al.l [2003k iLewis et al.l 
12006b . In particular, our ability to determine the precise value 
of the primordial spectral index of scalar fluctuations, n s , one 
of the key paramete rs to learn more about inflation (e.g. see 
iKomatsu et all |2010|. for recent constraints), could be severely 
compromised if the detailed physi cs of cosmological r ecom- 
bination are not understood (e.g. iRubino-Martm et ID HoTE 
for recent discussion). Over the past few year, this fact has 
motivated a lar ge number of works on the physics of recombina- 




Dubrovich & Grachev! 120051 : IChluba & Sunvacv 



Kholupenko & Ivanchik! |2006|; 



Wong & Scott 



2007 



Hiratal 



20081; IC 
t al.l |20C 



Switzer & Hirata 
Rubino-Martm et alj |2008|; 



Karshenboim & Ivanov 120081 : 
2008 1: Ijentschural 12009] ; lLabzowskv et 

20101) , all with the aim to get ready for the analysis of CMB data 



Chluba & Sunyaev 
2009; lGrin& Hirata 



from Planck, Act, Spt and in the future from Cmbpol. 

In this paper we want to discuss the effect of recombinations 
to highly excited states (n > 100) in hydrogen during cosmologi- 
cal recombination (z ~ 1000) on both the cosmological recombi- 
nation history and the recombination spectrum. In particular, we 
want to focus on the effects associated with the detailed evolu- 
tion of the populations in the angular momentum sub-states of hy- 
drogen, a process that already has received so me attention earlier 
jRubifio-Martm et al.ll2006r . lChluba et al.ll2007l) . However, the pre- 
vious computations were limited to m odels with 100 /-resolved 
shells in hydrogen dChluba etai]|2007t) . amounting to a total of 
5050 levels. As pointed out there, in order to obtain converged re- 
sults for the free electron fraction at low redshifts one has to include 
~ 200 - 300 shells, a task that requires refined numerical methods 
and significant improvements of the earlier recombination code. 

Recently, iGrin & Hiratl feOld) strongly advanced such com- 
putations, including up to 250 shells, or a total of ~ 3 1 000 sepa- 
rate states, during hydrogen recombination. They conclude, that for 
the computations of the CMB temperature and polarization power 
spectra the results for the cosmological recombination history are 
converged at the ~ 0.5cr at Fisher-matrix level when including 128 
/-resolved shells in the atomic model of hydrogen. Here we also 
reach a similar conclusion, using a completely different numerical 
method and independent recombination code, showing that for the 
Planck data analysis ~ 100 shells should already be sufficient at 
large multipoles /. Nevertheless, due to huge improvements in the 
performance of our recombination code, it is now possible to solve 
a recombination history for 100 shells in about 11 min on 8 cores, 
and even 350-shell computations only take a little longer than a day. 
Therefore, it will be easy to account for this correction providing 



1 For conditions in our Universe, the free electron and proton number den- 
sities never freeze out completely, but at low redshifts they enter a period 



a new training set for the multi-dimensional regression code Rico 
dFendt et alj|2009h . 

We also present res ults for the cosmological recombination 
spect ru m from hydrogen 1 Dubrovichl 19751; Dubrovich & Stolvarovl 



19951: IKholupenko et all 



2005 



Rubino-Martm et all l2006t 



Chluba & Sunvaevl l2006al ; ISunvaev & Chlubal l2009h " and provid 



more detailed computations which include the effect of collisional 
processes. In particular, we develop a new solver for the coupled 
system of ordinary differential equations (ODEs), which can be 
easily adapted to other problems, e.g. for computations of nuclear 
networks appearing in supernova and star form ation calculations 
(e.g. see iTimmesI U999I : Iffix & Meveil l200ot) . or for chemical 



Shapiro & Kang 


Il987t lAnninos et alj|l997l; iTegmark et al. 


1 19971; 


Abeletal. 1997 


; Gnedin et al.l 20091, and references therein) and 



l2008h . 

The paper is structured as follows: in Sect. [2] we give a few 
details and references about the cosmological recombination prob- 
lem. For more bas ic o verview we refer th e inter ested reader to 
ISeager et al.l d2000h and lSunvaev & Chlubd J2009h . In Sect. [3] we 
provide some details about the new recombination code and the 
ODE solver that we developed. This section is rather technical, and 
only meant for the interested reader. In Sect.f4]we discuss the results 
for the cosmological recombination history and the recombination 
spectrum. There we also briefly discuss the effects of collisions dur- 
ing recombination, and conclude in Sect. [5] 



2 DESCRIPTION OF THE RECOMBINATION PROBLEM 

Although the current version of our recombination code 
dChluba & Sunvaevl |2010|) allows us to include several additional 
processes that already have been shown to be important to pre- 
cise computations of the re combination history dFendt et al.|[2009l: 
iRubino-Martfn et al] l201Ct) . here we want to focus on the ef- 
fect of recombinations to highly excited level (n > 100) in hy- 
drogen, taking into account the detailed evolution of the pop- 
ulations in the angular momentum levels in each shell with 
principle quantum number n. We therefore restrict ourselves to 
a minimal model of hydrogen and helium, that does not in- 
clude any of t he detailed radiative transfer effects, such as pho- 
ton feedback dChluba & Sunvaevl 120071; ISwitzer & Hiratal 120081 : 



IKholupenko et all |2010|: IChluba & Sunvaevl 1201017 or Lym an-tt 



over which they evolve very slowly. 



diffusion dHirata & Forbes! 120091 : IChluba & Sunvaevl l2009tj) . De- 
tails about our atomic mo del f or hydrogen c an be found in 
IRubino-Martfn et alj d2006l) a nd IChluba et al.l d2007|). Th e he- 
lium model is explaine d in IRubino-Martfn et ah 20081) and 
IChluba & Sunvaevl d2010h. Details about the setup of the rate equa- 
tions can be found in lSeager et al. (2000). 

However, for the computations carried out in this paper, we do 
include the acceleration of helium recomb ination caused by the ab- 
sorption of resonant photons by hydrog en ( Kholupenk o et al.l2007l : 
ISwitzer & Hiratdl2008l ; IRubino-Martfn et al.ll2008h . This process 
allows helium to finish recombining until redshift z ~ 1700, so that 
accounting for this process is very important for the initial condi- 
tion of the hydrogen recombination problem, which at z ~ 1650 can 
be approximated using the Saha-equations. In most of our compu- 
tations we use this setup, however, when computing the cosmolog- 
ical recombination spectrum we evolve both hydrogen and helium 
starting at z = 3400 with S aha- values for their populations. 
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2.1 Numerical computations of the recombination and 
photoionization rates 

One of the crucial and also rather time-consuming aspects of the re- 
combination problem is the separate computation of the photoion- 
ization and recombination rates for each considered level. Here 
several points are important. For high levels the /-dependence of 
the photoionization cross section, <r, c (v) for level i = (n, Z), is 
very strong. This implies that for large n both recombination and 
photoionization are mainly occurring through the lo w-/ states, as 
the G aunt-factors drop very strongly with / (e.g. see lChluba et al.l 
120071) . Also the energy-dependence of the photoionization cross 
section becomes exponential for states with n » 1 and / ~ n, so 
that care must be taken when integrating the cross sections over the 
CMB blackbody spectrum. 

In addition, for high n stated the effect of stimulated recom- 
binations caused by the presence of CMB blackbod y photons has 
to be included. As shown earlier dChluba et al]|2007t) . not only the 
recombination rate for each level increases due to stimulated re- 
combinations, but also the /-dependence of the recombination co- 
efficients within a highly excited shell is strongly affected. 

Furthermore, for large n we encountered difficulties with the 
recursion relations for th e photoionization cross section given by 
IStorev & Hummed i ll99ll) . when going far away (several hundred 
times) from the threshold frequency v, c . In this case the initial con- 
dition for the recursion relation numerically became zero, so that 
the cross section is effectively zero for all / states. We avoid this 
problem by rescaling the terms in the recursion formulae with the 
n th root of the initial condition, which one could easily compute. 
With this the computations of the cross sections becomes stable up 
to very large n (~ 1000) and /, and very large distances from v, c . 

In order to choose the range of integration most efficiently, we 
identified the frequency v s > v, c at which cr, c v 2 becomes extremely 
smal0 (~ 10~ 30 - 10" 40 of the value at the threshold). For high n 
and / states, this occurred over only a few threshold energies, while 
for the low / states we found a much wider range. However, we 
typically limited the frequency interval to v, c < v < 10 8 v, c , and 
in the Wien tail of the CMB blackbody spectrum we in addition 
utilized the exponential cutoff for hv » kT y . 

Also, to avoid the time-consuming recursive computation of 
the cross sections at every time-step, we tabulated them on a dense 
grid and then use spline interpolation. This procedure accelerates 
the computations of the recombination and photoionization rates 
by a large factor, without significant loss of accuracy. 

To carry out the recombination integrals we implemented a 
ful ly adaptiv e integ rator based on the integration formulae given 
by IPattersor] dl968l) . This method is based on Gaussian quadra- 
ture rules, but in contrast to simpler Gaussian rules (e.g. Gauss- 
Kronrad) these formulae are fully nested, so that function eval- 
uations of every previous subdivision can be reused. We also at- 
tempted a scheme bas ed on Chebychev integration rules, but even- 
tually the formulae bv lPattersonl {l968) performed better. We com- 

2 It is straightforward to show that in our Universe stimulated recombi- 
nations become ~ 10 times stronger than spontaneous recombinations for 
levels with principal quantum number n > 23 [ yjgg] 

3 Since the integrand of the photonionization coefficient, fjfc, in the CMB 
blackbody field scales as <x, c v 2 /[e' ,v /*' r r — 1] < (kT-y/h)<Tj c v oc v 2 , while 
<r, c v 2 oc v , the latter provides a rather conservative estimate for the be- 
haviour of the integrand with frequency. Note that only the ratio to the value 
at the threshold frequency is important, so that cr, c v 2 drops slower towards 
larger v than <r, c v. 



pared our results for the recombination and photoionization rates 
with those obtained using NAG-integrators, and found excellent 
agreement (to the level of precision which we used for the inte- 
gration; this was normally ~ 10~ 8 - 10~ 9 in relative terms). 



3 THE NEW MULTI-LEVEL RECOMBINATION CODE 

Based on our previous multi-level h ydrogen and helium recombina- 
tion code dChluba & Sunvaevl2010h . we developed a new solver for 
the system of coupled ordinary differential equations. In particular, 
we replaced all routines that were used from the commercial Nag 
librar)Q so that our new code now is completely non-commercial, 
and parts of it could be parallelized. 

For the development of this new solver on a single processor 
machine three points turned out to be very important: (i) because 
of the vastly different timescales involved in the evolution of the 
populations of the excited states in hydrogen and helium, we use a 
stiffly stable algorithms with adaptive step-size control; (ii) because 
of the size of the Jacobian matrix for the equation system it is im- 
portant to make use of its sparseness for both, storage and matrix 
operations; (iii) it is also important to use analytic expressions for 
the Jacobian matrix when possible, in order to achieve sufficient 
numerical stability and accuracy, and also to make the code faster. 

After managing these aspects we also parallelized parts of our 
code using OpenMP. However the largest boost in our performance 
(a factor ~ 150 in comparison to our previous version) is obtained 
because of the improvements in connection with (ii) and (iii). Be- 
low we provide some more details on each point. 

3.1 The stiffly stable ODE solver 

The system of rate equations describing the cosmological recom- 
bination problem is very stiff. Even in the effective three-level 
approximation, wh ich is used in the implementation of Recfast 
dSeageretalJll999l) . because of the vastly different timescales on 
which the electron temperature and the population of the hydrogen 
ground state adjust their values, the stiffness of the equations is al- 
ready so large that a stiffly stable solver is required to achieve high 
accuracy and good performance. 

Below we explain the essential parts of the ODE solver we 
developed for the recombination problem. This solver was imple- 
mented in a general way, so that it also can be applied to other prob- 
lems, as already pointed out in the introduction. We would also like 
to mention, that in comparison to the commercial Nag sparse ODE 
solver routine D02NJF we were able to achieve a slightly better 
performance. This is likely due to the higher order of the implicit 
equations for the time-stepping (6th order instead of 5th), the ef- 
ficient use of sparseness in the equation setup, and the direct and 
complete ability to tune all solver parameters independently. Fur- 
thermore, with the Nag solver we were unable to solve the recom- 
bination problem for more than 150 shells. 

3.1.1 Gear's method with adaptive time-step 

For stiff problems implicit methods are known to be numerically 
more desirable than explici t methods, becaus e of both performance 
and stability issues (e.g. see lPress et al .1 1992b . We chose an implicit 

4 See http://www.nag.co.uk/numeric/ 
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Gear's method lGeaj[T97ll) to evolve the solution to the recombi- 
nation problem from one time /, to the next r j+1 . This method is 
related to the backward diff erence formulae (BDF) methods (see 
ICurtiss & Hirschfeldej|l952h . 

The differential equation system for the recombination prob- 
lem can be cast into the form 
dy 
At 

where y = y(t) is the (high-dimensional) solution vector of the 
problem at time t. Denoting the solution at time r, as y, the implicit 
equation for the time-step using the Gear's method reads 

5 



f(t,yl 



(l) 



y M = hAy M +^K k y t - k , 

k=0 



(2) 



with y M = dy, + 1 / At = f(t M ,y i+1 ), h = t M - t u and y,_ k being the 
solution to the recombination problem at previous times £;_t < 

It was shown by Gear i n 1971 that only up to 6th order such 
equation is stiffly stable (see lGearlll97ll) , so that we truncated the 
sum in Eq. ((2) at k = 5. The coefficients A and k, can be determined 
using the Taylor-expansion of the solution y i+ i. Since initially only 
the solution at t is known, and because we wanted to allow both 
changing order and variable time-step we re-derived these coeffi- 
cients up to sixth order. They can be found i n Appendix |Al For 
constant time-step they resemble those given in lAnitid feOOaT 



3.1.2 Solving the non-linear equation system 

Equation {2} defines a system of non-linear equations that deter- 
mines the solution yi+i at time r, + i. One can linearize this equation 
assuming that one has a guess for the solution y? +1 which is c7o.se to 
the correct one. Inserting yj' +1 on the right hand side of Eq. {2} one 
obtains >f +[ . This then leads to 



a-hf}J f ]6y'; +l =y'; +l -yl r (3) 

where Jf is the Jacobian matrix of the system Eq. |QJ evaluated at 
y = yf + ! , and Sy p + , = yf + Y -yf + { , where y p *l is the solution to Eq. ©. 
Note that in general 6y p +l j= y p +l - y? v After obtaining <5yf , from 
Eq. ([3} one can use yj^ 1 = y p +l + Sy p j+l as new initial condition 
for the problem and then iterate until convergence is reached. We 
usually use a combination of relative and absolute error control to 
check for convergence. Especially for the free electron fraction, the 
electron temperature, and the number of electrons in the ground 
states of hydrogen and helium, tight settings (e ~ 1(T 8 ) for the 
relative accuracies were necessary. 

To solve for 6y p +l it is not useful to explicitly invert the matrix 
A = [1 - hpjf]. Even though 1 - hfijf might be very sparse, 
the inverse [1 - hp J/] -1 likely will be dense. Also, inversion will 
only lead to numerically meaningful results for well-conditioned 
matrices. In addition, for the recombination problem the sparseness 
is of order percent (see Sect. |3,2| >, so that simple multiplications b = 
A x are fast, while one expects a very large number of operations to 
invert [1-AjSJ/]. 

We therefore choose an iterative schem e, based on a stabi - 
lized bi-conjugate gradient method (e.g. see iBarrett et al. IE991). 
This method generates two orthogonal sequences of vectors for the 
matrix A. These vector sequences are the residuals of the iterations, 
which are iterated until convergence is reached. As preconditioning 
we simply use the inverse diagonal elements of the Jacobian. For 
the recombination problem this is sufficient. 

In order to avoid many expensive evaluations of the Jacobian 
matrix, we furthermore make use of Broyden's method dBrovdenl 



1 19651 : lGadFl979l : IPress et al.ll 19921) . Here we assume that in each 
update only the non-zero matrix elements should be changed, so 
that the sparseness of the matrix is preserved during the whole run. 
We find that this approximation works very well for the recombi- 
nation problem. 

To obtain an initial guess for the solution at each time-step we 
use a simple extrapolation method based on the solution at previous 
time-steps. This allows us to write 



k=0 



(4) 



where the coefficients y, are also given in Appendix lAl 



3.1.3 Ordering of the equations 

We tried several orderings of the equations, however, none of them 
lead to significant improvements of the performance. In particu- 
lar, we also tried the ordering of the hydrogen states suggested by 
lGrin&Hirat3 ( l201(]|) . without any apparent benefit. In their order- 
ing, states with the same angular quantum number / are grouped 
and then ordered in the sequence / = 0, 1, 2, ... n mllx - 1, where n max 
denotes the largest shell that was included. In our final implemen- 
tation we use 



y(0 



T, 
X 1 
X< 



HI 
Hoi 



(5) 



where the populations of the hydrogen states, X ( H1 , are ordered like 
Is, 2s, 2p, 3s, 3p, 3d, 4s, 4p, etc. and the helium states, X ; HeI , are or- 
dered in the same manner, first for the singlet and subsequently for 
the triplet atom. This is t he same ordering used in earl i er implemen- 
tation s of the problem teubino-Martfn et all 1 20061 : IChluba et"aTI 
120071) . The distribution of elements in the Jacobian matrix for or- 
dering A (as in lRubino-Marti'n et alj200rj : IChluba et al.l2007l) . and 
ordering B (like in I Grin & Hir ata 2010) are illustrated in Fig.[T] 



3.1.4 Quasi-stationary approximation for the excited states 

The excited states (n > 1) adjust their populations on time-scales 
much shorter than the expansion time. It is therefore possible to 
simplify the differential equation system using the quasi-stationary 
approximation for levels with n > 1. We tried such approach, but 
found that in that case the performance of our solver was reduced. 

We also tried assuming that only for levels with n » 1 the 
quasi-stationary approximation was applicable, again with no gain 
in performance. We therefore used the full system of ODE's in 
our computations. However, we would like to note that even in the 
quasi-stationary approximation for levels with n > 1 we could ob- 
tain accurate recombination histories and also recombination spec- 
tra. For the recombination problem the quasi-stationary approxima- 
tion seem valid for the excited states, however, the performance of 
our ODE solver was not suggesting such an approach. 

Roughly speaking, the quasi-stationary approximation 
amounts simply to dropping the 1 term in Eq. (3} relative to the 
hfj Jf and re-normalizing with respect to h. When considered in 
this way, it is clear that the majority of the effort already results 
from solving a large system of non-linear equations. Furthermore, 
leaving the 1 term amounts to a physical regularization should the 
Jacobian become formally singular. 
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Figure 1. Structure of the Jacobian matrix for different cases: top panel - 
ordering A; middle panel - ordering B; lower panel - ordering B with 1- 
changing collisions. In all cases 20 shells for hydrogen and 5 shells for he- 
lium were included. The equations related to helium mainly lead to entries 
in the lower right corner of the Jacobian. The sparseness is about 10%. 



3.2 Sparseness of the Jacobian matrix and its setup 

The size of the Jacobian matrix for the ODE system scales like 
~ n 4 /4 with the total number of shells, n, that are included into the 
atomic model of hydrogen. For example, the recombination model 
for a 100-shell hydrogen atom with resolved angular momentum 
quantum numbers leads to a system of ~ 5050 differential equa- 
tions, and hence a Jocobian with ~ 25 million entries. For a 350- 
shell hydrogen model one already deals with about ~ 61 000 ODEs 
and the full Jacobian has ~ 3.7 billion entries. 

However, because of the dipole selection rules, the transition 
matrix is rather sparse. Simple estimates show that the number of 
non-zero elements scales ~ j« 3 , so that the matrix is sparse at a 
level of percent for 100 shells^. One can find examples for the struc- 
ture of the Jacobian in Fig.[Tj 

Because of the scaling of the number of non-zero elements 
with n the Jacobian becomes even more sparse for larger n. In com- 
putations, it therefore is important to make use of this sparseness, 
in both storage and matrix operations. For this purpose, we adapted 
the SparseLib++ Library^ to store the Jacobian matrix. For the 
purpose of parallelization we used the compressed-column format. 
Also, the required linear algebra operations were implemented us- 
ing this library. Here only operations with non-zero elements were 
performed, so that the efficiency is very large as compared to full 
matrix routines. 

To compute the Jacobian of the system Eq. (0, we use both 
analytic and numerical derivatives. For all the bound-bound dipole 
transitions of hydrogen, we use fully analytic expressions, while 
for helium at this stage we always use numerical derivatives. These 
are computed with a two-point central difference formula, which 
is second-order accurate in the choser0 5X. Also we compute the 
derivative of the recombination rates with respect to T e analytically, 
leading to an additional recombination integral in the matrix setup. 
A typical evaluation of the Jacobian for ~ 100 shells in our current 
implementation requires about ~ 6 seconds on a standard single- 
processor machine. With parallelization using OpenMP we were 
able to gain a factor of 3 - 4 for this evaluation on 8 cores. 

We also comment, that we usually solved for the free electron 
fraction using X e = 1 - HX, HI + / He i - Z^, H , instead of the dif- 
ferential equation. This is possible due to particle conservation. In 
our current solver we simply replaced the ODE for the free electron 
fraction by the above algebraic equation. 

3.3 Performance of the code 

The performance of our new recombination code initially is mainly 
limited by the speed of the matrix vector operation Ax and the 
setup of the recombination and photoionization rates. 

To reduce the effort in connection with the recombination inte- 
grals, w e use the tabulation scheme described in lChluba & Sunvaevl 
d2009bl) . We again confirm the precision of this procedure, by 
comparing with computations in which all the recombination 
rates were explicitly evaluated at each time-step. Also we tab- 
ulate the changes in the escape probabilities of the Hei 2'Pi - 
1 1 So and 2 3 P[ - l'Sp resonances caused by the presence of neu- 
tral hydrogen dKholupenko et alj |2007| ; ISwitzer & Hiratd 120081: 

5 For ~ 270 shells the sparseness drops below one percent. 

6 See http://math.nist. gov/sparselib++/ 

7 We also tried higher order formulae, but the tradeoff caused by the addi- 
tional number of function evaluations was too large. With the second-order 
formula we normally used SX ~ 10~ 6 X. 
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Table 1. Performance of the recombination code. n raax denotes the included 
number of hydrogen shells, and « eq gives the number of hydrogen levels. 
'Run A' is executed on a standard single-processor machine (MacBook Pro, 
2.4 GHz Intel Core 2 Duo, 3 GB 667 MHz DDR2 SDRAM), while for 'Run 
A (S m)' and 'Run B (S m)' we used nodes of the SciNet GPC supercom- 
puter, each with two 2.53 GHz quad-core Intel Nehalem E5540s, and 16GB 
of 1066 MHz DDR3 SDRAM, where m gives the number of cores that 
were used. For additional details on the parameters for the different runs 
see Sect. 13. 31 



"max 




Run A 


Run A (S 1) 


Run A (S 8) 


Run B (S 8) 


100 


5050 


52 min 


33 min 


1 1 min 


25 min 


150 


11325 


3.4 h 


2.0 h 


1.3 h 


2.8 h 


200 


20100 


9.7 h 


5.5 h 


3.6 h 


8.0 h 


250 


31375 


22 h 


12 h 


8.0 h 


17 h 


300 


45150 




23 h 


14 h 


35 h 


350 


61425 




42 h 


27 h 


63 h 



Rubino-Martm et al using the results of a calculation for 

a 5-shell hydrogen with 5-shell helium model. With this setup, for 
a 100-shell-hydrogen atom about 40% of all time is spent for mul- 
tiplications A x, while about ~ 20% of time is spent computing the 
recombination integrals. For a larger number of shells the contribu- 
tion from the multiplications A x strongly increases (e.g. reaching 
~ 70% for 200 shells), so that parallelization of this part is ben- 
eficial. However, by parallelization we currently only achieve an 
additional factor of ~ 2 in comparison to a single-processor ma- 
chine. 

To measure the performance of our code we compared runs 
with different settings for the accuracy. If one is only interested in 
the cosmological recombination history, it in principle is possible 
to run the code in a faster mode. The helium recombination history 
can already be computed rather precisely with 5 shells for hydrogen 
and 5 shells for helium. It is then possible to start the computation 
of hydrogen recombination just after helium has become neutral, 
including a much larger amount of hydrogen levels. For the spec- 
tral distortions from hydrogen on the other hand it is important to 
evolve the whole system starting well before helium recombina- 
tion, since even hydrogen is emitting some amount of photons in 
every transition because of the reprocessing of helium photons. 

In TableQ]we show a few examples regarding the performance 
of our code. For Run A we use a simplified helium recombination 
history, and only start the full hydrogen recombination calculation 
at z = 1650 evolving everything until z = 200. For Run B we start 
at redshift z = 3400 and evolve both hydrogen and helium until 
z = 200. Due to memory restrictions we did not compute cases 
with more than 250 shells on a single-processor machine. 

As can be seen from Table [T] for Run A and n mlix > 200 on 
a single-processor machine the computational time scales roughly 
oc ;j^ 3 6 with the total number of hydrogen shells, or like oc n^ q 7_L8 
with the total number of equations. Starting the computation at red- 
shift z = 3400 instead of 1650 take about twice as long. On a single 
core our code seems to be a factor of ~ 10 fasted than the one of 
iGrin & Hiratal J2010h . 

On SciNet's GPC systenfl we find an immediate performance 
gain of a factor of approximately 1.6 on a single processor, despite 
the very small change in clock speed; this is almost certainly due 

8 Their run for 200 shells hydrogen starting at redshift z ~ 1606 takes about 
4 days on a single-processor machine.. 

9 See http://www.scinet.utoronto.ca/ 
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n - 125 




200 400 600 800 1000 1200 1400 

Figure 2. Correction to the ionization history for different number of hy- 
drogen shells as a function of redshift. Shown is the relative difference with 
respect to Recfast (see text for details). One can see that at high redshifts 
(z ^ 900) the dependence on n max is very small, while at low redshifts the 
results seem to start converging for 350 shells. 



to the improved memory bandwidth of the newer Nehalem proces- 
sors. The SparseLib++ matrix operations were parallelized with 
OpenMP; the matrix-vector multiply was decomposed by counting 
the number of matrix non-zeros and forcing a static decomposition 
which split the number of non-zeros approximately equally over 
the number of threads; other matrix operations were parallelized 
with dynamic partitioning over rows. In addition, we parallelized 
the Jacobian matrix setup, including the computations of the re- 
combination and photoionization integrals needed for updates of 
the interpolation tables. This allowed us to gain another factor of 
~ 2 using 8 cores, so that eventually a computation of the recombi- 
nation history and recombination spectrum was about 2-3 times 
faster than on a standard single-processor machine. For Run B on 
8 cores and for « raax > 200 the computational time again scaled 
roughly oc /ij lax 3 5 with the total number of hydrogen shells, or like 
oc n' q 7 ~ 18 with the total number of equations. 

We would like to mention that one could further speed the 
computations of the recombination history up when ignoring accu- 
racy for the cosmological recombination spectrum. With this one 
would likely be able to gain another factor ~ 2 - 3. For example, 
we always limited our step-size to Az < 3, to obtain a resolved 
representation of the recombination spectrum. This could easily be 
increased twice. Another possibility is to allow overshooting be- 
yond the desired redshift point, and then interpolating the solution 
where it is requested. This will diminish the accuracy for the spec- 
trum, but still leads to precise results for the recombination history. 
However, at this point we did not follow this idea any further. 



4 RESULTS 

In this section we present our results for the cosmological recom- 
bination history and the cosmological recombination spectrum. We 
first discuss the effect of the completeness of the atomic model and 
then present some first results in connections with collisional pro- 
cesses, which we incorporate using simple approximations. 
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Figure 3. Correction to the ionization history at different redshifts as a func- 
tion of the the total number of hydrogen shells. 




Figure 4. Total recombination coefficient, a tat = 2 n> ; O'nl, as a function of 
ft max and for different redshifts. For simplicity we assumed that T c = T y , 
and we also included the effect of stimulated recombinations in the ambient 
CMB blackbody field. 



4.1 The cosmological recombination history 

In this section we present the results for the corrections to the 
cosmological rec ombination history . We compared with the Rec- 
fast vl.4.2 code JWong et"al]|2008l) . but excluding the corrections 
to the helium recombin ation history and g etting rid of the switches 
in the ODE system (see lFendt et alj2009l for details). We used the 
standard hydrogen fudge factor, / H = 1.14. In Fig. [2] we present 
the results for the modification of the free electron fraction as a 
function of redshift including a different number of hydrogen shells 
(see Sect. 14.1.21 for a more detailed discussion on the shape of the 
correction). One can clearly see that at high redshifts, close to the 
maximum of the Thomson visibility function (z > 1000), the de- 
pendence of the correction on the number of hydrogen shells is 
already rather small. On the other hand at low redshifts, it is im- 
portant to include shells up to rc max ~ 300 - 350 to obtain fully 
converged results for the ionization history. 

This can be seen even more clearly in Fig. [3] which shows the 
dependence of the correction at fixed redshifts z on n raax . Appar- 
ently, at z ~ 200 the change to the recombination history starts to 
converge for n raax ~ 300 - 350, while at z > 900 already n max ~ 100 
is sufficient. As mentioned in the introduction, this is because at 




Figure 5. Changes to the CMB temperature and polarization power spectra 
when going from 100 shells to 350 shells. T he curves were obtained using 
a modified version of Cmbeasy tooranl2005l) . 



z <; 800 - 900 the dynamics of recombination is strongly controlled 
by the escape of photons from the Lyman-o- resonance and the net 
two-photon decay-rate of the 2s level, so that increasing the ef- 
fective recombination rate cannot affec t the ionization his tory very 
much. This was already pointed out in iFendt et all d2009h , m con- 
nection with the effect of the Recfast hydrogen fudge factor on the 
recombination history. 

On the other hand, at low redshift, the completeness of the 
atomic model is still important, and the 'bottle-neck' for recombi- 
nation is not only set by the 2s- Is transition rate or the Lyman-o- 
channel, but also by the capture rate of electrons from the contin- 
uum. However, at the level of percent to the correction at z ~ 200 
the atomic model seems to become complete for n max ~ 350. 

We would like to note that the simple total recombination co- 
efficient, tUtoi = Yinj a ni, continues to grow for n max > 350 (see 
Fig-EJ- Although a tBt does not capture any of the dynamical effects 
(e.g. net radiative transition rates and escape probabilities; possible 
collisional ionizations for very high levels) that define the effective 
recombination coefficient, Fig. [4] shows that in principle it is pos- 
sible to increase the effective recombination rate, when the bottle- 
necks of recombination set by the effective rate of transitions to the 
ground state will be modified. Such kind of modification could for 
example be achieved by collisions that mix different / sub-states, as 
we will discuss in Sect. 14.31 

We would also like to mention that until now we have not 
checked the convergence of the results at lower redshifts (z < 200). 
There the chemistry of the Universe will also become important. 
In addition, large uncertainties in the reionization physics will af- 
fect our ability to carry out precise computations. Resolving these 
questions is beyond the scope of this paper. 



4.1.1 Effect on the CMB power spectra 

For the analysis of Planck data, only the changes in the CMB tem- 
perature and polarization power spectra really matter. In Fig. [5] we 
show the correction that is obtained when going from 100 to 350 
shells. This correction is very small, only reaching AC//C; ~ 0.05% 
at / = 3000, and also remaining well below the cosmic vari- 
ance limit at small /'s. Therefore, one does not expect any im- 
portant changes to the biases in n s and fib/! 2 obtained earlier by 
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Rubino-Martm et al. I i boiol) using the results for a 110-shell hy- 
drogen models. Regarding the Planck data analysis ~ 100 shells 
already seem to be sufficient. Nevertheless, due to the large im- 
provement in the performance of our recombination code, it will 
be rather easy to account for the full correctio n, e.g. providing an 
updated training set for Rico dFendt et alj|2009h 

Although the final correction to the free electron fraction at 
z ~ 200 decreases by a factor ~ 1.8 when going from 100 to 
350 shells (see Fig. |2j, the modification in the CMB power spec- 
tra at large multipoles / is much smaller, simply because the main 
effect is connected with a change in the total optical depth r, 
while the position of the maximum (z ~ 1100) and the width 
(Az ~ 200) of the visibility function remain practically unchanged. 
However, the value of t could still be affected at the level of per- 
cent by this modification. One therefore expects a change in the 
EE power spectrum at low multipoles /, close to the region that 
is also affected by reionizat ion (e.g. see iHaiman & Holderl [20031 ; 
IColombo & Pierpaolil 120091 and references therein). We indeed 
found a AC//C/ ~ -0.8% correction to the EE power spectrum 
at low multipoles (/ < 10), while the TT power spectrum remained 
nearly unaltered. However, our computations did not include any 
detailed reionization model, so that it is difficult to estimate the 
impact of this correction regarding the value of t. Similarly, one 
will have to include possibl e effects of dark matter annihilation 



20091; ISlatver et al.l 120091: ICirelli et al.l 


120091; fautsi et alj 20091: 


Kanzaki et al.l201Cl 


) or decaying particles dChen & Kamionkowskil 


200% Zhang et al. 2007) into such considerations. This is bevond 



the scope of this paper. 



4.1.2 Possible approaches for improvements to Recfast 

The shape of the correction to the free electron fraction for our 350 
shell computation (cf. Fig.[2]and also the solid/black line in Fig.[6jl, 
suggests that in comparison with Recfast there are (at least) two 
different regimes: the region at z > 800 and the low redshift region 
at z < 800. At z > 800 the correction grows with a larger slope than 
at z < 800. Also the correction shows some small local maximum 
at z ~ 720 and a local minimum at z ~ 590. What are the likely 
physical reasons for this behaviour? 

The are several aspects of the recombination problem that cur- 
rently cannot be captured by Recfast. First one might think about 
the effective recombination coefficient, a B , which in Recfast is 
probably slightly overesti mated. There, a B is mo delled using the 
fitting formulae given by IPequignot et al.l J 1 99 ll> pl us some ad- 
dition al overall hydrogen fudge factor / H = 1.14 dSeager et al .1 
1 19991) . By adjusting the hydrogen fudge factor one can compen- 
sate part of the co rrection, as already pointed out earlier (e.g. see 
lFendtetal]|2009l) . Here we found that for the 350 shell results 
/h = 1.110-1.115 captures the amplitude and rough scaling of 
the modification around z ~ 900, but overestimates the effect at 
z < 800 by a factor of 1.5 - 2.0 (see Fig. [6] blue/dashed curve). 
Also, by simply changing /h it is impossible to reproduce the local 
maximum or minimum at z ~ 720 and z ~ 590. 

There are two additional aspects of the recombination problem 
that change at z ~ 800: the temperature of the matter starts to depart 
significantly from the photon temperature. This will modify details 



10 Here the subscript 'B' stand for case B recombination , where direct re- 
combinations to the ground state are excluded (e.g. see iBaker & Menzel 
ll938l ; lHummeJl994) . 



2h ;md 2j> modelled 
mJqK.'iiduLtlv + tuning 1 
2h ;md 2p modelled 
independently + tuning II 




Figure 6. Correction to the ionization history with respect to Recfast for 
different recombination models. The solid/black curve shows the result 
from our multi-level recombination code for « max = 350. The other curves 
were computed using Recfast, with different modifications relative to the 
standard case (see text for details). 



of the interplay between the photon field with the hydrogen atom, 
so that cfa becomes a function of both electron and photon tempera- 
tur e. This aspect cannot be captured with the simple formulae given 
bv lPequignot et all dl99lh . and without additional computations it 
is hard to estimate the importance of this effect. However, as we 
show below, the temperature scaling of a B indeed changes the low 
redshift behaviour of the correction. 

In addition, around the same redshift the 2s- Is two-photon 
channel again becomes less important than the Lyman-a transition, 
since the optical depth in the Lyman-a resonance drops. This has 
more delicate consequences, which currently cannot be captured 
by Recfast. Since in the derivation of the Recfast equations it is 
assumed that the 2s and 2p state are in full statistical equilibrium 
{Ni p = 3 A?2s), every electron that reaches either the 2s or the 2p 
state is immediately shared with the other. It was already shown 
earlier, that this assumption is not valid at th e end of recombina - 
tion, since collisional processes are inefficient dChluba et al .120071) . 
Avoiding this approximation therefore would allow to distinguish 
between electrons that reach the 2s state and the 2p state. A simple 
derivation (using the quasi-stationary assumption for the evolution 
of the 2s and 2p populations) yields a slightly modified Recfast 
equation for hydrogen 

dX H 1 

^ = ^7JTT){^[^-^^] 

+ C 2v [x e N p a 2p -3X ls ^2 P ]}, (6) 

where H(z) is the Hubble expansion factor, X c and X [s are the free 
electron fraction and ls-population, N p = [1 - X ls ] A'h, and 7V H is 
the total number of hydrogen nuclei in the Universe. Furthermore, 
the recombination and photonionization rates for the 2s and 2p state 
are denoted as a, and f} h respectively. Due to the presence of CMB 
photons one also obtains the factor £ = exp(-/iV2i /k T y ), where 
v 2 i ~ 2.47 x 10 15 Hz is the Lyman-cf transition frequency, and the 
2s and 2p inhibition factors are given by 



C 2s 



A 2s i s +/3 2s 



(7a) 



(7b) 



Here A 2s i s is the 2s- Is two-photon decay rate, and the effective 
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Lyman-o- transition rate is given by A 



2pls 



PsA 2 , 



where P$ is the Sobolev escape probability of the Lyman-cy reso- 
nance, and A^pis ~ 6.27 x 10 8 s is the Lyman-a transition rate. 

The important point about Eq. @ is that one has to spec- 
ify the partial recombination rates to the 2s and 2p state. Using 
detailed balance one then also obtains the partial photonionzation 
rates just like in Recfast. This allows to account for more details in 
the microphysics of the multi-level cascade. Given a B , one would 
normally assume a 2s s ;Sb and a 2s ~ | (Yb, so that fj 2s w ^/?b 
and y6 2p ~ \Pb, where Pb is the photoionization coefficient in the 
normal Recfast case. However, the detailed dynamics of recombi- 
nation connected with the cascade of electrons from higher levels 
to the 2s and 2p state (these are included into the computation of 
qtb), can render this approximation crude, in particular at the end of 
recombination. One aspect that is connected with this, is that the 2s 
and 2p state will departure from full statistical equilibr ium, as seen 
earlier dRubino-Martrn et al.ll2006l : IChluba et ai]|2007l) . 

To demonstrate the principle possibilities of Eq. we show 
the resulting correction with respect to the normal Recfast, choos- 
ing / H = 1.135, a 2s = 0.59 q-b, and a 2s = 0.41 a B (Fig. [6] curve 
label with 'tuning F). As can be seen, there now is a maximum at 
z ~ 800, which is actually produced by the fact that the Lyman-a 
transitions again starts controlling the recombination process. Note 
that in general the ratio y = a 2i /a 2p should be a function of red- 
shift, so that one can accommodate more details of by computing 
a 2s and a 2p using the results of our full recombination code. 

Finally, we also computed the correction for /h = 1.135, 
a 2s = 0.59 q-b, and a 2s = 0.41 q-b, but using &b = ctsiTy) instead 
of Ob = <*B(Xe) (Fig. [6] curve label with 'tuning IF). Recombi- 
nations physically are controlled by both the electron temperature 
and the photon temperature. The latter enters, because CMB photon 
strongly control the radiative cascade within the hydrogen atom, 
and we just wanted to show, that varying the temperature depen- 
dence of &b in fact has most effect at low redshifts. We therefore 
expect that a detailed computation of a 2s and a 2p from the full re- 
sults of our recombination code, may allow us to capture more as- 
pects of the microphysics. 

Along the same line, one may consider including more lev- 
els into the Recfast code (e.g. 3s, 3d, 3p). Again one could de- 
termine the effective recombination rates from the results of our 
recombination code for each included level. Such approach may al- 
low to capture in detail some of the feedback processes within th e 
Lyman-series dChluba & Sunvaevll2007l ; IKholupenko et alj|2010h . 
In addition, simple modifications should allow to include the ef- 
fect of stimulated 2s- Is two-photon dec ays and the reabsorption of 
Lyman-a photons in the ls-2s channel dChluba & Sunvaevj|2006bl ; 
IKholupenko & Ivanchikll2006i) . In combination with our new ODE 
solver, such extensions will probably not degrade the speed of Rec- 
fast very much, so that it can still be used in the analysis of future 
CMB data. We plan to investigate these possibilities in the future. 




Figure 7. The bound-bound cosmological recombination spectrum of hy- 
drogen. The lower panel shows a zoom in on the frequency range around 
v ~ 1 GHz. Free-free absorption and collisions were not included. 



tion from hydrogen still increases by a factor ~ 1.6 at v ~ 0.1 GHz, 
when going from 100 to 350 shells. However, at this level the spec- 
trum seems to converge. From the observational point of view the 
region v > 1 GHz is much more interesting. There it may become 
possible to measure the v ariable component of the recombination 
spectrum in the future (see lSunvaev & Chluball2009l for overview). 
At v > 1 GHz the recombination spectrum seems to converge at the 
level of a few percent and better, when including 350 shells. How- 
ever, this statement is only true regarding the completeness of the 
atomic model for hydrogen. The effect of /-changing collisions still 
alters the shape of the recombination spectrum at low frequencies 
dChluba et alj|2007t) . as we will further discuss in the next section. 
Also at frequen cies v < 0. 1 GHz f ree-free absorption should be- 
come important dChluba et al.l2007t) . so that accurate computations 
in this region are more involved. 



4.2 The cosmological recombination spectrum of hydrogen 
at low frequencies 

Including more than 100 shells in the computation of the hydrogen 
recombination spectrum should m ainly affect the low frequency 
part of the recombination spectrum dChluba et al .l2007h . This is be- 
cause including more shells enables more electrons to pass through 
highly excited states, emitting additional low frequency photons in 
the cascade towards lower levels. 

In Fig.|7]one can observe this effect. The recombination radia- 



4.3 The effect of collisions 

Already in our earlier work dChluba et al1l2007h we discussed the 
effect of collisions on the cosmological recombination spectrum 
and the ionization history of our Universe. However, because of 
computational restriction in that implementation we were unable to 
include more than 100 shells for hydrogen. However, it is expected 
that in a diluted plasma such as in our Universe in the epoch of re- 
combination, the effects of collisions will become more important 
for higher shells. The most important aspect of this problem is that 
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Figure 8. Ratio of the total /-changing collision rate caused by electron and 
proton impact to the total radiative transition rate in the CMB blackbody 
radiation field for different shells of hydrogen and at different redshifts as a 
function of /. Here / raax = n - 1 . The kink at low / is due to transitions out 
of the rep states, which depopulate fastest. 



collisions can lead to mixing of / and n-states in the Rydberg levels. 
One expects that /-mixing becomes important for n > «/_ m i x , and 
then 72-mixing and mixing with the continuum starts at n > M n _mi x , 
where n„_ m i x » "/-mix s> 1, since /7-changing collisions (which re- 
quire transfer of energy) are less effective than /-changing collision. 

Here we investigate the effect of collisions on both the ioniza- 
tion history and the cosmological recombination spectrum for up to 
300 shells. We include /-changing collisions, n-changing collisions, 



n = 100 
n = 200 
n = 300 
n = 400 
n = 500 



500 



1000 



1500 



2000 
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Figure 9. Ratio of the total collisional ionization rate induced by electron 
impact to the total photoionization rate caused by absorption of photons 
from the CMB blackbody for different shells of hydrogen. At z > 1000 
collisional ionization is always faster than photoionization for n > 400 - 
500. Note however, that collisional ionization is much slower (by a factor 
of ~ 10 2 - 10 ) than bound-bound radiative transitions until much larger 
redshifts (z > 3000 - 4000). Collisional ionization therefore never seems to 
be very important for the recombination problem. 



and collisional ionizations by electron, proton and a-particle (only 
important at z > 1 80 0) impact in our computations, as described in 
IChluba et all ( l2007h . In that work, the main colli sional rates were 
compu t ed using simple approxim a tions given by van Regemortej 
jl962h : |Pengellv & Seatonl dl964»:lBrockMiurstl Jl97ll) . For addT 
tional overview see also lMashonkinal jl99^X 

A more detailed treatment of collisional processes is far be- 
yond the scope of this paper, but given that the effects on the low 
frequency recombination spectrum are significant, it may become 
important to refine these computations. With our current estimates 
for the collisional rates, the dynamics of recombination do not seem 
to be affected at a level that is very important for the CMB power 
spectra, however, we would like to point out that the accuracy of the 
used collisional rates easily allows for factors of a few. Therefore, 
it will be very important to refine the current calculations of the 
collisional rates, in order to give a final answer for their relevance 
during recombination. 

Additionally, at low redshifts, collisions with neutral hydrogen 
atoms could start to become important (see iMihailov et al]l2004l . 
for recent computations of rate coefficients). However, we defer 
this problem to future work. 



4.3.1 Importance of the different collisional processes 

To understand the effect of collisions and at which n they are ex- 
pected to become important, it is illustrative to compare the main 
collisional rates with the radiative rates in the CMB blackbody ra- 
diation field. Here in particular stimulated emission and photon ex- 
citation are important, and we included both into our computations. 

In Fig. [8] we show the comparison of the /-changing collision 
rate with the total radiative rate, including emission, excitation and 
ionization processes, as a function of / for different n. We com- 
puted the rates for the mixture of electrons and protons as obtained 
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with the Recfast code. At the considered redshifts helium is al- 
ready completely neutral, and hence was excluded. As one can see, 
at z ~ 200 only the very high Rydberg states are expected to be 
mixed over / by collisions, while states with n < 200 should start to 
drop out of full statistical equilibriurrFI (SE). At z ~ 600 one finds 
"/-mix ~ 150, and at z ~ 1000 even states with n ~ 100 should start 
to be completely mixed. 

This actually suggests another simplification of our recombi- 
nation code. For levels with n > 300 - 400 it should always be 
possible to treat them as if they are in full SE. This should allow to 
add many more shells to the recombination problem, since it will 
be possible to obtain the solution for the populations within a given 
shell n > «/_ ra j x with only one additional differential equation per 
shell. However, at this stage we have not followed this possibility 
any further. 

We also checked the effect of collisional ionization. In Fig. [9] 
we show the importance of collisional ionization by electron im- 
pact relative to the photoionization rate as a function of redshift. 
Here it is important that according to the used approximations the 
collisional ionization rate is directly proportional to the photoion- 
ization cross sections. Therefore collisional ionization, like pho- 
toionization, is only effective from low-/ states (/// max ~ 0. 1 — 0.3). 
Even if collisional ionization can become comparable to the normal 
photoionization rate, in particular for the high /-states normal radia- 
tive bound-bound transitions will dictate the evolution of the levels. 
Collisional ionization therefore never seems to be very important 
for the recombination problem. Our numerical computations also 
confirm this statement, where only the effect of /-changing colli- 
sions seems to be notable. 

Regarding collisional excitations and de-excitations, we found 
that for tt-transi tions (An = 1) among e xcited levels, with the for- 
mulae given by Ivan Regemorteil i 19621) at z ~ 1100 already for 
n ~ 30 - 50 these rates become comparable to the radiative tran- 
sitions. However, here it is important that for most of the levels 
in some given shell n the a-transitions are not defining the popu- 
lations of the levels, so that the effect on the recombination pro- 
cess remains small up to much larger n. Therefore, also collisional 
bound-bound transitions should be of minor importance for the re- 
combination problem, but one will have to refine the modelling of 
these rates for detailed predictions of the low-frequency recombi- 
nation spectrum (v < 1 GHz). There energy changing collisional 
transitions and collisional ionizations are expected to suppress the 
emission of photons during the recombination epoch. We find such 
behaviour when artificially enhancing the collisional excitation and 
ionization rates. 

We would like to note, that collisional excitation and de- 
excitation, and collisional ionization at some very high n are ex- 
pected to push the populations of levels towards an equilibrium 
with the continuum. Since collisions are mediated by particle im- 
pact, the temperature of the electrons will be important for these 
levels. However, in our computations up to /7 raax = 300 the pop- 
ulations of the levels were always more close to Boltzmann equi- 
librium with the lower states, where the CMB blackbody temper- 
ature is important. Nevertheless, even for those shells that started 
to be /-mixed in our computations the deviations from Boltzmann 
equilibrium with lower states were still rather large. All this sug- 
gests that for additional simplification of the recombination prob- 
lem, the probably most straightforward approach will be to repre- 



In full SE one has N„i = (21 + l)N ns for the level populations within a 
given shells n. 
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Figure 10. Correction to the ionization history caused by collisional pro- 
cesses for different number of hydrogen shells as a function of redshift. 

sent /-mixed shells with only one additional differential equations, 
but to account for the evolution with n in detail. 

4.3.2 Effect on the cosmological ionization history 

In Fig.[l0]we present the correction to the ionization history caused 
by collisions. According to our computations /-changing collisions 
are dominating the corrections. At high redshifts, collisions do not 
affect the recombination history much. There the effective recom- 
bination rate is already saturated, so that mixing among the high / 
states, which should lead to an increase in the effective recombina- 
tion coefficient, does not result in any modification. At low redshifts 
(z s; 1000), the increase in the effective recombination coefficient 
mediated by /-changing collisions, leads to a small acceleration of 
recombination. 

At this point we did not treat the effect of collisions for more 
than 300 shells. However, one does expect some additional changes 
when going to a larger number of shells, as the fully mixed levels 
will continue to accelerate recombination. But, here it will be very 
important to refine the calculations of the collisional rates, as with 
the current accuracy they may easily be up to a few times off. We 
therefore, stopped at this point and will return to this problem, once 
we have better estimates for the collisional rates. 

4.3.3 Effect on the cosmological recombination spectrum 

In Fig.QT]we show the effect of collisions on the cosmological re- 
combination spectrum. In particular at low frequencies (v < 1 GHz) 
collisions modify the recombination radiation in comparison to the 
case without collisions (see Fig. [7}. Due to /-changing collisions, 
more electrons reach the high n and high / states. From these levels 
only transitions with small An are possible, so that more photons 
are emitted at low frequencies, explaining the enhancement there. 
We would like to mention, that for the considered cases, n-changing 
collisions and collisional ionizations are still not very important. 

We also note that at this stage the uncertainty in the collisional 
rates does not allow to compute precise templates for the recombi- 
nation spectrum at v < 1 GHz. Also in connection with the ioniza- 
tion history, the final modification caused by collisions cannot be 
given at this point. It will therefore be important to refine the treat- 
ment of collisional processes. We plan to investigate this problem 
in the near future. 
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Figure 11. Effect of collisions on the bound-bound cosmological recom- 
bination spectrum of hydrogen. For the upper panel /-changing collisions, 
collisional excitation and ionization were included. In the lower panel we 
show the direct comparison of the obtained spectra for n max = 300 with and 
without collisions included. Free-free absorption was not included. 



5 CONCLUSIONS 

In this paper we compute the cosmological recombination history 
including a total of ~ 61 000 states into the hydrogen atom. For 
this purpose we developed a new ODE solver that allows us to treat 
large systems of very stiff differential equations. This solver was 
implemented in a general way, so that it also can be applied to 
other problems, e.g. in computations of nuclear networks appear- 
ing in supernova and star formation calculations, or for chemical 
networks important during reionization and the dark ages. 

We discussed the effect of recombinations to the high Rydberg 
states in hydrogen on the dynamics of recombination and the re- 
combination spectrum. Our computations indicate that for the cos- 
mological ionization history at redshifts z > 200 the atomic model 
for hydrogen becomes complete for n max ~ 300 - 350. At z ~ 200 
we find a modification of AN S /N e ~ 1.6% when including 350 
shells. This is about a factor of 1.8 smaller than the result for a 
100-shell computations (see Sect. 14. 11 and Fig.[2]for more details). 
Nevertheless, the effect on the CMB temperature and polarization 
power spectra at large multipoles is small when going from 100 
to 350 shells (see Fig. [5J. From the point of view of Planck data 
analysis ~ 100 s hells seem to be already eno ugh, and the conclu- 
sions reached by iRubino-Martfn et all d201fj|) for the biases to n s 
and Cl b h 2 caused by the detailed physics of recombination should 
remain practically the same. However, the total value of t could 
still be affected at a significant level by the discussed process. 



Due to the huge improvements in the performance of our re- 
combination code it is possible to compute the recombination his- 
tory for 350 shells in about a day (see Table[T|l. Although according 
to our computations 100 shells seem to be sufficient for Planck data 
analysis, we still plan to provide an updated training set for Rico 
dFendt et alj|2009h , which then could be used for precise computa- 
tions of the CMB power spectra. Also for the final code compari- 
son, this type of updated training set may be useful. Alternatively, 
we plan to investigate possibilities for a refined treatment within a 
Recfast type of scheme, as mentioned in Sect. 14. 1.21 

The inclusion of more shells to hydrogen also leads to an in- 
crease of the amplitude of the cosmological recombination spec- 
trum at low frequencies, since recombinations to highly excited 
states allow additional electrons to emit photons in transitions with 
small An. The main results for the recombination spectrum are pre- 
sented in Sect. |4.2| and in particular Fig. [7] At frequencies around 
and above ~ 1 GHz, the recombination spectrum converges at the 
level of percent or better, when including 350 shells. Such preci- 
sion will become important in the future when using measurements 
of the cosmologica l recombination spectrum to constrain cosmo- 
logical parameters dSunvaev & Chlubal2009h . Nevertheless, addi- 
tional non-standard processes may als o be important here, e.g. du e 
to pre-recombinational energy release JChluba & Sunvaevll2009bh . 
or annihilation of dark matter particles dChluball201fj|) . 

We also discussed the effect of collision on both the dynam- 
ics of recombination and the recombination spectrum (Sect. 14.3.21 
and I4.3.3t . Our analysis suggests that at z > 200 levels with 
n > 300 - 400 will always be completely mixed over /. Colli- 
sional ionizations and excitations, however, seem to become im- 
portant only for much larger n. If states are completely mixed over 
/ it will, in principle, become possible to add many more shells to 
the recombination problem, since only one additional equation per 
shell will be required above ni_ m j x . We plan to investigate this possi- 
bility in the near future, however, at this point it will be more urgent 
to refine the computations of collisional rate coefficients. Here we 
only used very rough approximations, common for computations 
in stellar astrophysics. However, both the cosmological recombi- 
nation spectrum and the recombination dynamics require updated 
and more precise computations of these rates. We hope that this 
work will motivate some experts in atomic physics to attack this 
complicated problem in the near future. 

With our current estimates we find a correction to the free elec- 
tron number density of AA^/A^ 8.8 x 10~ 4 at z ~ 700, which 

is mainly caused by /-changing collisions with protons. However, 
this result could be off by factors of a few because of the uncer- 
tainties in the used collisional rates. Also, for the final answer one 
will probably have to include more shells into the computation, as 
the convergence of this correction is very slow with n max . Our cur- 
rent computations with collisions were limited to n max = 300, how- 
ever, according to our estimates n max ~ 300 - 400 will likely be 
necessary. For computations of the CMB power spectra it will be 
important to get this correction right. 
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APPENDIX A: COEFFICIENTS FOR THE GEAR'S 
METHOD WITH VARIABLE ORDER AND TIME-STEP 

To obtain the coefficients for the variable time-step implicit Gear's 
formula Eq. ® and the extrapolation formula Eq. © it is conve- 
nient to define 



Af;-* = tf- k - t, 

Atj-k 
Af; +] ' 



Pk 



(Al) 
(A2) 



Here k = 0, 1, 2, 3, 4 when the solution at 4 previous and the current 
time tj are available. For k = one has Af, = and p k = 0. 



Al Coefficients for implicit Gear's formulae 

Using the definitions given above one finds 
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One has to calculate the k, starting with k 5 . To use a lower order 
Gear's formula / < 5 one has to set k, = for i > I. Note that in our 
notation n*[-] = 1 f° r ^ ^ /• 



A2 Coefficients for extrapolation 

Using the definitions given above one finds 



yo 



Yj 



1 - 
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Pi 



5 

k=l 

i-iyf\^-tp^t\ P — 
UPi-pt fitt L\Pj-p* 



for j < 5. One has to calculate the y j starting with y 5 . To calculate 
the extrapolation based on 1 < / < 5 previous time-steps one has to 
set yj for j > I equal to zero. Note that again nit—] = 1 f° r k > j. 
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